In [1]:
import ps_util as ps
import numpy as np
import matplotlib.pyplot as plt
import plotly
plotly.offline.init_notebook_mode() 
%matplotlib inline

Introduction to Photometric Stereo

Read in the Beethoven Data, you will be working with primarily 3 data artifacts:

3D Array: Images (m, n, k)

2D binary mask : (m,n)

array S : light vectors

In [2]:
#read in the file which has I,M,S where I is the images, M is the mask, and S is the light vectors.
V=ps.read_data_file('Beethoven')

Step One: Get aquainted with the data

Images, mask, light vectors

In [7]:
# Take the I,M,S and store them into variables
Iv=V[0]
mv= V[1]
sv= V[2]

plt.figure(figsize=(7,7))

plt.imshow(Iv)
plt.show()

plt.figure(figsize=(7,7))
plt.imshow(mv)
plt.show()

plt.figure(figsize=(7,7))
plt.imshow(sv)
plt.show()
In [13]:
#displaying the images
for i in range (0,3):
    plt.figure(figsize=(7,7))
    plt.imshow(Iv[:,:,i])
    plt.show()

Step 2: Applying Mask to Images

In [14]:
# use the mask_image to calculate the relevant values within the mask (where values > 0) 
nz=np.where(mv)
len(nz[0])
Out[14]:
28898
In [15]:
mv.shape
Out[15]:
(256, 256)
In [16]:
# Take each individual image and apply the mask to it, then reshape the outupt into a (3,nz) array. 
a=Iv[:,:,0]
b=Iv[:,:,1]
c=Iv[:,:,2]

ai=a[nz]
bi=b[nz]
ci=c[nz]

J=np.row_stack((ai,bi,ci))
J.shape
Out[16]:
(3, 28898)

Step 3: Applying Lamberts Equation

I is the intensity of each pixel in the image, p is the albedo of each pixel, N is the surface normal and S is the lighting source with the assumption that it is constant over each image.

Take the dot product of the Image times the inverse of the light vector to obtain the Modulated Normal Field.

In [35]:
#take the light intesnity vector, take its inverse and store it in a variable to use in Lamberts Equation
# M = S^-1 * J 
SV=np.linalg.inv(sv)
SV.shape
print(SV) 
[[ 1.62948467 -3.51636848  1.70279567]
 [ 1.95828261 -0.13769666 -1.75540327]
 [ 0.36274532  0.32154769  0.31388367]]

Step 4: Calculate the Albedo Modulated Normalized Field

In [36]:
#Take the dot product of the two matrices in order to solve for M 
Mb=np.dot(SV,J)
In [37]:
Mb.shape
Out[37]:
(3, 28898)

Step 5: Calculate the Albedo

Albedo: material light absorption property

In [38]:
#Normalize the M to obtain the albedo
albedo=np.linalg.norm(Mb, axis=0)
albedo.shape
Out[38]:
(28898,)
In [39]:
#Store the albedo back into an empty numpy array and visualize the output. 
Beto=np.zeros((256,256))
Beto[nz]=albedo
fig, ax = plt.subplots(figsize=(7,7))
plt.imshow(Beto)
plt.title("Albedo for Beethoven in 2D")
plt.axis('off')
Out[39]:
(-0.5, 255.5, 255.5, -0.5)
In [40]:
#Albedo with a different color scheme
plt.figure(figsize=(7,7))
plt.imshow(Beto, cmap='binary')
Out[40]:
<matplotlib.image.AxesImage at 0x10ef03a90>

Calculate the normal fields

In [41]:
#Calculate the normal fields n1, n2, n3
n1=np.divide(Mb[0],albedo)
n2=np.divide(Mb[1],albedo)
n3=np.divide(Mb[2],albedo)
In [42]:
Nbus=np.vstack((n1,n2,n3))
In [43]:
#Take the normals and store them into an empty numpy
be1=np.zeros((256,256))
be2=np.zeros((256,256))
be3=np.ones((256,256))
In [44]:
be1[nz]=n1
be2[nz]=n2
be3[nz]=n3

Betov=np.dstack((be1,be2,be3))

Betov.shape
Out[44]:
(256, 256, 3)
In [45]:
#code to use to plot the image through plotly
x = np.where(Betov)[0]
y = np.where(Betov)[1]

x.shape
y.shape
Out[45]:
(123332,)

Display Normal Field

In [46]:
#display images 
plt.figure(figsize=(6,6))
plt.imshow(be1)
plt.title("Normal Field 1")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be2)
plt.title("Normal Field 2")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be3)
plt.title("Normal Field 3")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(Betov)
plt.title("All Normals")
plt.show()
In [47]:
#display images 
plt.figure(figsize=(6,6))
plt.imshow(be1, cmap='binary')
plt.title("Normal Field 1")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be2, cmap= 'binary')
plt.title("Normal Field 2")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be3, cmap='binary')
plt.title("Normal Field 3")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(Betov, cmap='binary')
plt.title("All Normals")
plt.show()

Step 8: Retrieve 3D

In [48]:
#calculate the z to plot using unbiased integrate 
z0=ps.unbiased_integrate(be1,be2,be3,mv)
z0.shape
Out[48]:
(256, 256)
In [49]:
ps.display_depth_matplotlib(z0)
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning:

elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison

In [50]:
#calculate the z to plot using simchony 
z1=ps.simchony_integrate(be1,be2,be3,mv)
In [81]:
ps.display_depth_matplotlib(z1)
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:584: RuntimeWarning: invalid value encountered in less
  cbook._putmask(xa, xa < 0.0, -1)
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1545: RuntimeWarning: invalid value encountered in greater
  intensity > 0),
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1550: RuntimeWarning: invalid value encountered in greater
  hsv[:, :, 2] = np.where(intensity > 0,
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1556: RuntimeWarning: invalid value encountered in less
  intensity < 0),
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1560: RuntimeWarning: invalid value encountered in less
  hsv[:, :, 2] = np.where(intensity < 0,
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):
In [87]:
#since matplotlib and mayavi were subpar options, we opted to graph using a plotly interactive graph
In [88]:
import plotly.plotly as py
from plotly.graph_objs import *
In [163]:
trace1 = Surface(
    z=z1,  # link the fxy 2d numpy array
    x=y,  # link 1d numpy array of x coords
    y=y,   # link 1d numpy array of y coords
    colorscale="Greys"
)
In [164]:
data = Data([trace1])
In [169]:
# Dictionary of style options for all axes
axis = dict(
    showbackground=True, # (!) show axis background
    backgroundcolor="rgb(204, 204, 204)", # set background color to grey
    gridcolor="rgb(255, 255, 255)",       # set grid line color
    zerolinecolor="rgb(255, 255, 255)",   # set zero grid line color
)

# Make a layout object
layout = Layout(
    title='Beethoven', # set plot title
    scene=Scene(  # (!) axes are part of a 'scene' in 3d plots
        xaxis=XAxis(axis, showticklabels=False), # set x-axis style
        yaxis=YAxis(axis, showticklabels=False), # set y-axis style
        zaxis=ZAxis(axis, showticklabels=False)  # set z-axis style
    )
)
In [170]:
# Make a figure object
fig = Figure(data=data, layout=layout)
In [171]:
#the graph is interactive and you can zoom in and out of the 3D figure
plotly.offline.iplot(fig)
Drawing...

Assignment:

Albedo Modulated Normal Field

Albedo

Normal Field

Depth -Z